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ABSTRACT 

In this work we have developed a new stochastic model for the fluctuations in lightcurves of accreting 
black holes. The model is based on a linear combination of stochastic processes and is also the solution 
to the linear diffusion equation perturbed by a spatially correlated noise field. This allows flexible 
modeling of the power spectral density (PSD), and we derive the likelihood function for the process, 
enabling one to estimate the parameters of the process, including break frequencies in the PSD. 
Our statistical technique is computationally efficient, unbiased by aliasing and red noise leak, and 
fully accounts for irregular sampling and measurement errors. We show that our stochastic model 
provides a good approximation to the X-ray lightcurves of galactic black holes, and the optical and 
X-ray lightcurves of AGN. We use the estimated time scales of our stochastic model to recover the 
correlation between characteristic time scale of the high frequency X-ray fluctuations and black hole 
mass for AGN, including two new 'detections' of the time scale for Fairall 9 and NGC 5548. We find a 
tight anti-correlation between the black hole mass and the amplitude of the driving noise field, which is 
proportional to the amplitude of the high frequency X-ray PSD, and we estimate that this parameter 
gives black hole mass estimates to within ^0.2 dex precision, potentially the most accurate method 
for AGN yet. We also find evidence that w 13% of AGN optical PSDs fall off flatter than l// 2 , and, 
similar to previous work, find that the optical fluctuations are more suppressed on short time scales 
compared to the X-rays, but are larger on long time scales, suggesting the optical fluctuations are not 
solely due to reprocessing of X-rays. 

Subject headings: accretion, accretion disks — galaxies: active — methods: data analysis — quasars: 
general 



1. INTRODUCTION 

Both galactic black holes (GBH) and active galac- 
tic nuclei (AGN), being black hole accretion systems, 
share a number of similarities in their spectroscopic 
and variability features despite being separated by sev- 
eral decades in mass. The major ingredients of the 
two systems appear to be an optically thick, geomet- 
rically thin accretion disk emitting thermal radiation, 
and an optically thin hot 'corona' emitting X-ray emis- 
sion with a power-law spectrum. Because the tem- 
perature of the disk scales with the black hole mass 

as kT oc Mgg , the thermal disk emits in the soft 
X-rays for GBHs, and in the optical/UV region for 
AGN. The source and geometry of the corona emis- 
sion is not known, but likely involves in verse compton 
emiss i on from a hot electron pla sma (e.g.. iShapiro et all 
119761 : lHaardt fc Maraschil 1 1 9 9 lh with a possible syn- 
chrot ron contribution from a jet {Markoff e t al.l I2001L 
12005ft . Both GBHs and AGN follow a correlation be- 
tween the blac k hole mass, radio luminosity, and X- 
ray luminosity dMerloni et al.ll2qol iKording et al.ll2006b 
lYuan et al.1l2009HGultekin et al.ll2009bft . leading some to 
suggest that the systems are largely self-similar, at least 
in low accretion r ate states, differin g only in mass and 
environment fe.g.. iFalcke et aLll2004f) . 

However, the comparison between GBHs and AGNs 
is not straightforward, as GBHs are known to ex- 
ist in different 'states', with the same source seen 
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cycling between the various states (for a review see 
iRemillard fc McClintock1l2006t lBeIIonj|2010ft . The spec- 
tral and timing properties of GBHs arc different in the 
different states. The most commonly observed states are 
the hard, soft, and very high state (VHS). In the hard 
state, the spectrum is dominated by the power-law hard 
X-ray emission, while in the soft state the X-ray flux 
increases and the spectrum contains a significant ther- 
mal disk component. In the VHS, the X-ray flux is very 
high and the spectrum is intermediate in hardness be- 
tween the soft and hard states. The power spectral den- 
sity (PSD, P(f)) of the X-ray fluctuations can be rea- 
sonably approximated by bending power-laws, P(f) oc 
1 / f a . For the soft and intermediate states, the power-law 
component contributes to the majority of the variabil- 
ity ( e.g., iDone fc Gierlifiskil 120051: iSobolewska fc Zvckil 
12006ft . while for the hard state the disk blackbody com- 
ponent dominates the soft X-ray variability, at least on 
time scales > 1 sec (jWilkinson fc Uttlevll2l)0flh . The PSD 
at high frequencies has a logarithmic slope a > 2, flat- 
tening to a slope of a ~ 1 below some break frequency. 
For GBH in the hard state and VHS, there is an addi- 
tional flattening to a = below a second lower break 
frequency, with the break frequencies being higher in the 
VHS. The accretion rate increases as one moves from the 
hard state to the VHS. 

For AGN the situation is less clear, due to the fact 
that the time scales involved for the state transitions 
are expected to scale upwards with black hole mass. 
As such, AGN, with their supermassive black holes 
(SMBH), have not been observed to unabmiguously 
undergo state transitions. Many AGN display SEDs 
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more characteristic of soft state GBH, in that there 
is a strong thermal component from t he di sk in the 
optical/UV region. iSobolewska eTail ((2009) suggest, 
based on comparisons of the AGN and GBH SEDs, that 
most AGN in current surveys are in the VHS (see also 
ISobole wska. Siemigino wska. fc Gierlinskyll201ll ). Almost 
all of the ~ 10 AGN with high quality X-ray lightcurves 
exhibit PS Ds similar to those seen in Cyg X-l in the sof t 
state fe.g.. lMarkowitz et al J 120031; IMcHardv et~aTl l2006). 
The only exception i s Akn 564, which exhibits a second 
low-f requency break ([Arevalo et a l. 2006; McH ardv et al l 
2007f) . Th e high accretion rate ( rh ~ 1. IRomano et al.l 
20041 ) led IMcHardv et~aTI ([20071 ) to suggest that this 



source is analagous to the VHS. However, the aver- 
age 2-10 keV photon i ndex of these sources is r < 2 
(jPapadakis et al.l 120091; ISobolewska fc Papadakii 120091) . 
which resembles the value observed for GBHs in the 
hard state; thus, there is some descrepancy between 
the spectral and timing classifications. The break fre- 
quencies for both GBHs and AGN are observed to anti- 
correlate with black hole mass (lUttlev fc McHardvil200l 
IMcHardv et al.|[2006f ). further strengthening the similar- 
ity between GBHs and AGN. The anti-correlation be- 
tween the break frequency and black hole mass has typ- 
ically been interpreted as being driven by a correlation 
between the size of the X-ray emitting region and black 
hole mass. A number of studies have also found an 
anti-correlation bet ween the amplitude of X-ray variabil- 
ity and Mbh (e . g. lLu~fcYul [20011 IBian fc Zhaol 12001 
Papadakisl l200l iNikolaiuk et al.l 120041; IQ'Neill et al l 
20051; iGierlinski et all 120081; iZhou et al.ll2010D . The ab- 



sence of AGN analogous to the hard state suggests that 
selection effects may be at work, as all GBHs with accre- 
tion rates less than a few per cent of Eddington are hard 
state objects (e.g., IMaccarong 12003?) . while most AGN 
surveys have only been able detect o bjects radiating at 
L/LEdd ^ 0-1 in large numbers (e.g ., iVestergaardl l2004t 
iTrump et al1l200l IKellv et al.ll20Tol) . 

While studies of AGN variability have many com- 
plications, they do have the advantage that the disk 
and corona emission can be cleanly separated. In con- 
trast, for GBH in the soft state the seperation is more 
difficult as both components emit in X-rays, and the 
disk com ponent does not exh ibit significant variabil- 
ity (e.g., IChurazov et al.l l200l"l ). Therefore, studies of 
the disk emission variability are more easily carried 
out for AGN. Previous studies of AGN optical vari- 
ability using well-sampled lightcurves have found that 
the optical variat i ons have dispersions of ~ 10-20% 
(|Kellv et all [200l iMacLeod etall l2010al ). with the X- 
rays varying more on the shorter time scales (e.g., 
' Ulrich et alj|1997t ICzernv et al.ll2003t ISmith fc Vaughanl 
20071: lArevalo et all 12001)7 and sometimes also on 
the longer time scales ([Breedt et"al"1 120101 ). In ad- 
dition, the optical PSD can be well described by a 
power law form P(f) oc I If 2 (IGiveon et all Il99l 
ICollier fc Peterson! [200lt ICzernv et al.n 20031 or, when 
the lightcurve i s long enough, by a Lorentzian cen- 
tered at zero (ICzernv et all 119991; IKellv et al.l 120091; 
IKozlowski et all 1201ft MacLeod et al.l I2010aj) ; Le.. the 
PSD is P(f) oc l// 2 , flattening to a constant below 
some break frequency. Similar to the X-ray PSD, the 
break frequency of the optical PSD increases with black 



hole mass (ICollier fc Peterson! |200H IKellv et all l200l 
IMacLeod et al.ll2010aD . although not as steeply. 

While variability studies of AGN and GBHs have been 
important for understanding the similarities between 
these two classes of objects, variability studies are im- 
portant for other reasons as well. For one, variability 
is one of the only observational tools available for prob- 
ing the disk viscosity, as the time scales of variability 
are expected to depend on v iscosit y, while the S EP is 
not (iSiemiginowska fc Czernvl[T98llStarling et al.ll200l 
iFrank. King, fc Raind 120021 ) . Understanding the disk 
viscosity is important for understanding the transfer and 
removal of angular momentum in the disk, which is 
fundamental to an understanding of the accretion flow. 
Variability is also a potentially important observational 
tool for constraining the geometry of the corona, with, 
for example, studies of quasi-periodic oscillations be- 
ing consistent with the disk evaporatin g into a hot in- 
ner flow at lower accretion rates (e .g., ICui et ail 119991 ; 
iRossi et all 12004 iDone et al.l 120071 ) . and the origin of 
the variability in the QPOs being in the hot corona 
in the soft state and the disk in the hard state (e.g., 
ISobolewska fc Z vcki 20 061). This is also c onsistent with 
the finding of iWilkinson fc Uttlevl (|2009D who showed 
that the disk is responsible for much of the variabil- 
ity in the hard state; however, there is still much work 
to be done, as it has not been conclusively shown that 
the source of variability originates in the same region 
as the emission. In addition, variability may be the 
most effective observational discriminator between the 
different accretion states, and thus may provide evi- 
dence of radiatively inefficient accretion flows (RIAFs) 
in AGN, which are believed to be associated with the 
hard state, should a hard state exist for AGN. A number 
of RIAF candidates exis t (e.g., E3E999; Qu ataert et all 
H999tlTrump et aDl2009T ) . but their spectral/timing state 
is unknown. Studies and confirmation of AGNs in the 
hard state are particularly important, as the har d state 
is associated with jet production in GBHs (e .g., iFenderl 
[200lt iFender et alTl200l iKording et all [20061 ) . and me- 
chanical feedback from these jets plays an important role 
in heating intracl uster gas and regulating the growth of 
massive galaxies (Croton et al.ll2006t iBower et al.l 120061 
ISiiacki et alll2007l) . 

Time series exhibiting PSDs of a 1// type are known 
as 'long-memory' processes, and there is a n extensive lit - 
erature on them (e.g., a good reference is iPalmal 120071 ). 
Almost all previous studies of variability in GBHs and 
AGN have relied on non-parameteric techniques, such 
as the periodogram or the structure function. However, 
there are a number of known difficulties in estimating 
the PSD or structure function non-pa rameterically (e.g., 
see lVaughan et al.l [2003: Pcs sahll2007T ). For one, the em- 
pirical estimate of the PSD, known as the periodogram, 
suffers from windowing effects due to the finite sampling 
of the lightcurve. These windowing effects include red 
noise leak and aliasing, which are caused by power leak- 
ing into the periodogram from time scales longer and 
shorter than the maximum and minimum time scales 
probed by the lightcurve, respectively. Red noise leak 
and aliasing distort the periodogram, making it poten- 
tially difficult to relate the observed periodogram to the 
true intrinsic PSD. Moreover, irregular sampling further 
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distorts the periodogram, although this distortion can 
be alleviat ed through the use of the Lomb-Scargle pe- 
riodogram (lLomblll976t lScargldll982t iHorne fc Baliunasl 
119861 (Zechmeister Kursterl 120091: iVio et al.1 120101 ) . The 
structure function is not immune to these problems and 
can si milarly be distorted, making its interpretation dif- 
ficult (E mmanoulopoulos et al.l 12010). These distortion 
problems can be significant for AGN especially, as their 
observed flux is often much fainter than GBHs and their 
lightcurves tend to be more irregular and sparsely sam- 
pled. 

M otivated by t hese p roblems, lUttlev et al.l (|2002t see 
also iDone et all (|1992f )) developed a powerful Monte 
Carlo technique for estimating the underlying PSD which 
accounts for the distorting effects of finite, and possibly 
irregular, sampling. The basic idea behind the technique 
is to use Monte Carlo simulations to calculate the ex- 
pected value of the periodogram as a function of the 
true underlying PSD, and then to find the true PSD 
which minimizes a y 2 goodness of fit measure between 
the observed periodogram and the expected one. Con- 
struction of confidence regi ons may be obtained t hroug h 
the procedure outlined by iMueller fc Madeiskil ([2009). 
This technique has the advantage that it may be used 
to fit any arbitrary PSD. However, it also has two dis- 
advantages. First, it is computationally intensive, espe- 
cially when constructing confidence regions. Second, the 
X 2 goodness of fit statistic used to fit the PSDs ignores 
the covariance in the periodogram among the frequency 
bins, and is not proportional to the log-likelihood of the 
perio dogram. As s uch, minimization of the x 2 statis- 
tic of lUttlev et al.l (|2002[ ). while effective, is unlikely to 
be the most efficient means of constraining the under- 
lying PSD. An alternative Bayesian approach based on 
an approximation to the likelihoo d function of th e pe- 
riodogram has been developed by IVaughanl (|2010D . An 
ad ditional likelihood-b ased approach has been developed 
by iMiller et al.l (|2010l) with the difference being that the 
likelihood function is calculated in the time domain. Un- 
fortunately these two alternative approaches do not com- 
pletely incorporate the distortion in the PSD due to fi- 
nite sampling of the lightcurve, such as that caused by 
re d noise leak, althou gh in principle the PSD model used 
bv lMiller et all (f2010h can be modified to correct for this. 
Indeed, it is very difficult to derive an analytic expres- 
sion for the likelihood function of the periodogram for 
a finite and irregularly sampled lightcurve, and in many 
cases may be impossible. 

An alternative and complementary approach to 
Fourier-based techniques is to model the lightcurve as 
a parameterized stochastic process, with the parameters 
of the model being related to the underlying PSD (or 
rather, the PSD being a function of the parameters of 
the process). Under this approach, the parameters of the 
model are estimated directly from the lightcurve itself; no 
Fourier transforms are performed, and thus there is no 
spectral distortion. Moreover, modeling the lightcurve 
in the time domain potentially can provide further in- 
sight on features in the power spectrum, such as break 
frequencies, as it is not always apparent how to inter- 
pr et the PSD. Time-d omain modeling was advoc ated for 
bv iKellv et al.l (l200l herea fter KBS09. see also iScargld 
(|1981D and lRvbicki fe Press! (|1992ft ) within the context of 
estimating the characteristic time scale of AGN optical 



variability (or equivalently, the break frequency of the 
optical PSD), and they developed a Bayesian approach 
for estimating the parameters. KBS09 modeled AGN op- 
tical lightcurves as Gaussian Ornstein-Uhlenbeck (OU) 
process, the power spectrum of which is a Lorentzian 
centered at zero, and showed that this process is consis- 
tent with the optical lightcurves of AGN for their sam- 
ple. Subsequent work has confirmed that the OU process 
provides a good description of AGN optical variability, 
at lea st on time scales between a few days and se veral 
years (jKozlowski et al.ll2010tlMacLeod et al.ll2010al) . and 
of blazar sub-mm variability (<Strom et al.ll2010D . More- 
over, the OU process provides a framework in which to 
measure the time lags betwe en the AGN br oad line re- 
gion and optical continuum dZu et al.ll2010f) and select 
quasars (IKozlowski et al.l 120101: IButler fc Blooml [20101: 
IMacLeod et al.H2010bD . ~ 

In this work we extend the method of KBS09 to enable 
more flexible modeling of the PSD of accreting black hole 
systems. We model the lightcurves as a linear expansion 
of OU processes, which enable us to approximate many 
of the features seen in the PSDs of GBHs and AGN, par- 
ticularly the bending power-law forms. We derive the 
likelihood function for this statistical model and perform 
statistical inference within a Bayesian framework, thus 
enabling one to calculate the probability distribution of 
the parameters, such as the break frequencies, given the 
data. Because our method is based on the likelihood 
function, it uses all of the information in the data. Fur- 
thermore, it fully accounts for measurement errors, ir- 
regularly sampling, red noise leak, and aliasing. Fitting 
is performed on the entire lightcurve simultaneously re- 
gardless of the sampling, thus enabling one to easily com- 
bine lightcurves obtained from different instruments and 
sampling time scales. These properties make our tech- 
nique particularly attractive for obtaining constraints on 
PSD features for poorly sampled lightcurves of faint ob- 
jects, as our Bayesian method uses all of the informa- 
tion in the data. Calculation of the likelihood function 
is computationally efficient, and we are able to calcu- 
late a maximum-likelihood estimate of the PSD in under 
a minutcQ. The primary disadvantage of our method is 
that, while flexible, it cannot fi t arbitrary PSDs; fo r this, 
one should use the method of lUttlev et al.l ([2002), or a 
combination of the two techniques. 

While the OU process, and similarly the mixture of 
OU processes, is useful for fitting PSDs, and thus quan- 
tifying variability, it is not always clear how to interpret 
the best -fit OU process astrop hysically, or the PSD in 
general. iTitarchuk et al.l (|2007D studied the linear diffu- 
sion equation for an accretion disk with a driving noise 
term as a model for the PSDs of accreting black holes, 
showing that features in the PSD depend on the viscosity 
of the accretion flow. In this work we will show how the 
mixture of OU processes arises as the solution to the lin- 
ear diffusion equation, perturbed by a random spatially- 
correlated white noise field. This interpretation of the 
mixed OU process may be considered to be among the 
'perturbation' class of astrophysical models for variabil- 
ity of GBHs and AGN, with the variability arising from 

1 The calculation was done on a lightcurve with ~ 3000 data 
points in IDL on a Mac Pro with two 3.2 GHz Quad-Core Intel 
Xeon processors 
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small random accretion rate perturbations that propa- 
gate inwards through the accretion flow, making the ob- 
served variations the prod uct of perturbati o ns that oc- 
curr a t larger radii (e.g . . iLvubarskiil 119971: iKing et all 
20041: lArevalo fe UttlevT 120061: Uaniuk fc Czernvl l2007t 



Titarchuk et alJl2007t ) . The perturbation class of models 



for variability also explains the correlation between flux 
and absol ute RMS variability seen in bot h GBHs and 
AGN (e.g. lUttlev fc McHardvll200lh . which lUttlev etafl 
(2005) show is very well approximated as being due to a 
Gaussian process on the logarithmic scale. 

Finally, we note that the perturbation class of mod- 
els for variability is only applicable to the broad-band 
flickering noise seen in GBH and AGN lightcurves, and 
we do not address the more catastrophic, and seem- 
ingly not stochastic, changes that have b een observed 
in the lightcurves of GRS 1915+105 (e.g.. iBeUoni et al.1 
2000). These may occur due to accretion disk insta- 
bilities, such as radiation pressure instabilities (e.g. , 
Shakura fc Sunvaevl 119761: iLightman fc Eardlevl Il974t 
Czernv et alJl2009D. or the thermal- viscous ionization in- 
stability (ILin fc Shielddll986l : ISiemiginowska et al.lll996t 
Uaniuk et al.l 12004). Furthermore, we stress that the as- 
trophysical interpretation of the mixed OU process as 
a model for the fluctuations from GBHs and AGN is 
based on the linear diffusion equation, which is surely an 
oversimplification, and MHD simulations must be per- 
formed for studying more physically motivated models 
for variability of accre t ing b la ck holes, as done b y, e.g 



Armitaee fc Reynolds! (12003ft iSchnittman et all (2006 
Moscibrodzka et all (120071) . iRevnolds fc Millerl (|2009D 
and lNoble fc Krolik] (|2009D . Rather, we study the mixed 



OU process as a solution to the stochastic linear diffusion 
equation in order to provide a guide for interpreting the 
best-fit model and PSD features of lightcurves, both real 
and simulated, as analytical solutions may be obtained. 

The format of the paper is as follows: In § [2] we de- 
scribe the mixed OU process statistical model for the 
lightcurves of GBHs and AGN. In §[3] we show that the 
mixed OU process is generically the solution to the lin- 
ear diffusion equation. In § 14.11 we derive the likelihood 
function of the mixed OU process parameters, and pos- 
terior probability distribution of the parameters given an 
observed lightcurve. In §[5] we apply the mixed OU pro- 
cess to an X-ray lightcurve of Cygnus X-l in the low/hard 
state, the X-ray lightcurves of 10 local well-studied AGN, 
and the optical lightcurves of a sample of AGN. Finally, 
in §|6]we summarize our results. 

2. THE STATISTICAL MODEL 

2.1. The Ornstein- Uhlenbeck Processes 

In this section we give a brief overview of the properties 
of the OU process that are relevant for o ur wo rk. Fur- 
ther details can be found in Kelly et al.l ([2009, see also 
IGillespie! ([19961 ) and IGardinerl (|2004D ). who refer to this 
process as a first-order continuous autoregressive process. 
The OU process, X(t), is a simple stochastic process by 
which the quantity of interest (say, the logarithm of the 
luminosity or accretion rate) responds to an input noise 
process with an exponential decay to its mean. Math- 
ematically, the OU process is defined by the following 
stochastic differential equation: 

dX(t) = -u) (X(t) - fx)dt + <;dW(t), w ,<;>0. (1) 



The parameters of the process are the characteristic fre- 
quency, wo, the mean value of the process, fi, and the am- 
plitude of the driving noise process, The term ? 2 has 
units of RMS 2 sec -1 , and thus <j gives the rate at which 
variability power is injected into the stochastic process 
X(t). The term W(t) denotes a Wiener process (i.e., a 
Brownian motion), and its derivative is white noise. Al- 
though in this work we will focus on the special case when 
W{t) is a Wiener process, the results outlined here for 
the OU process, such as the form of its PSD, are valid for 
a more general class of stochastic processes called Levy 
processes. For completeness, we give a brief description 
of Levy processes in the appendix. It is apparent from 
setting <; = in Equation ([TJ that X(t) decays to its 
mean value with an e-folding time scale of r = \/lu$. 
The time scale r is often called the relaxation time scale 
of the process; in this work we will refer to r as the char- 
acteristic time scale of the process, as r is related to the 
break frequency in the power spectrum of X(t). 

When the driving noise, dW(t), is white, the OU pro- 
cess is stationary and Markov. A stationary process is 
one whose joint probability distribution does not change 
when shifted in time, and a Markov process is one whose 
future states only depends on the current state. In addi- 
tion, if W(t) has zero mean and unit variance the auto- 
covariance function of the OU process is 



? 2 T 

Rou{t) = — e 



\t\/r 



(2) 



It follows from setting t = in Equation ([2]) that the vari- 
ance of the OU process is T 2 t/2, and that the autocorre- 
lation function has an exponential decay with t = 1/loq 
being the decorrelation time scale. The power spectrum 
of the OU process is given by the Fourier transform of 
the autocovariance function: 



Pou{u) 
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277 UJq 



(3) 



The power spectrum of an OU process is flat for fre- 
quencies lj <C ujq, and decays as 1/lu 2 for frequencies 
uj 3> luq. Hence the association of r as a characteristic 
time scale of the process: the OU process, X(t), resem- 
bles white noise on time scales long compared to r, and 
resembles red noise on time scales short compared to r. 
Note that if one calculates the power spectrum in terms 
of ordinary frequencies instead of angular frequencies, 
then r = l/(27r/o), where /o = ljq/(2tt). This implies 
that simply fitting the break in the power spectrum as a 
function of ordinary frequency / will overestimate r by 
a factor of 2ir. 

2.2. Mixtures of Ornstein- Uhlenbeck Processes 

We can build on the OU process to develop more flexi- 
ble models that enable modeling of more complex power 
spectra. In this work we consider a mixture of OU pro- 
cesses with independent driving noises for constructing 
a process with a power spectrum of the form exhibited 
by the X-ray lightcurves of GBH and AGN. We define a 
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discrete mixture of OU processes as 



Y M {t)=H + 



M 

E« 



iXi(t). 



(4) 



Here, c%,...,cm is a set of mixing weights, and 
Xi(t), . . . , Xjvf(t) is a set of OU processes with character- 
istic frequencies ojx , . . . , ujm and driving noise amplitudes 
Ci , . . . , ?m • A continuous version of Equation (j4]) is de- 
scribed in the Appendix, although we do not use it in 
this work. The autocovariance function of Ym {i) is 



M 



C 2 Z 2 T- 
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(5) 



where we have used the fact that r,- = 1 /u)i . The sum of 
exponentials in Equation falls off more slowly than 
a single exponential function, and therefore the autocor- 
relation function of the mixed OU process falls off more 
slowly than the autocorrelation function of the individual 
OU processes. In other words, the mixed OU process ex- 
hibits longer range dependency than a single OU process, 
and therefore is better suited for modeling lightcurves 
with long time scale dependencies. 
The power spectrum of Ym (t) is 



P 



Y,M 



(*) = E 



1 



2tt 



uj- 



(6) 



From Equations ([5]) and ([6]) it is apparent that the con- 
tribution of an individual OU process to the variability 
amplitude of a lightcurve depends on the product qq, 
and thus the two parameters are degenerate. Therefore, 
for simplicity in this work we only consider the case where 
all of the individual OU processes have the same value 
of?. 

In Figure Q] we plot a mixed OU process sampled at 
the three different time intervals. The mixing weights 
for this example were chosen such that the power spec- 
trum of the simulated lightcurve was flat below a low- 
frequency break, cjl, decays as 1// above the low- 
frequency break, and then steepened to 1/f 2 above a 
high-frequency break, ujh, similar to what is seen in the 
X-ray lightcurves of accreting black hole systems. Thus, 
the simulated lightcurves probe the 'white' noise, 'pink' 
noise, and 'red' noise regions of the PSD. We constructed 
these lightcurves assuming that the logarithm of the flux 
follows a mixed OU process with /j, = 0, ? = 1, th = 0.1 = 
1/loh, and t£ = 10 = 1/u>l, and that the driving noise is 
Gaussian white noise. In general, we model the logarithm 
of the flux as a Gaussian mixed-OU process, as this is 
consistent with the flux-rms correlation seen in the X-ray 
variations of GBHs and AGN (jUttlev fc McHardvll200lh 
and the log-normal distribution of fluxes for Cygnus X-l 
(jUttlev et al.l 12005). In addition, the same seed for the 
random number simulator was used in order enable a 
more direct comparison between the lightcurves. As is 
apparent, the mixed OU process can appear very differ- 
ent depending on which time scales are probed by the 
observation. 

Equation (IB7[) in the Appendix demonstrates that the 
power spectrum of the continuous mixed OU process, 
with mixing function given by Equation (|B2|) . exhibits 
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Figure 1. Simulated lightcurves illustrating the mixed OU pro- 
cess. The lightcurves were generated using fi = 0, c = X,tjj = 
0.1, Ti = 10, independent Gaussian white driving noises, and a 
common random number generator seed. The three lightcurves 
probe the white noise (flat) part of the power spectrum (top), the 
'pink' noise (1//) part of the power spectrum (middle), and the red 
noise (l// 2 ) part of the power spectrum. The mixed OU process 
can look very different depending on the time scales probed by the 
observations. 



the same behavior that is seen in the power spectrum 
of the X-ray lightcurves of GBHs and AGN. This there- 
fore suggests that the mixed OU process well represents 
the fluctuations in the X-ray lightcurves of GBH and 
AGN, and may hold clues to the physical origin of the 
X-ray fluctuations. Moreover, we have identified a pro- 
cess which gives an explicit connection between the break 
frequencies in the X-ray power spectra and characteristic 
time scales. This is important, because it is not always 
apparent how to connect a feature in the power spec- 
trum with behavior in the time-domain. For example, 
we note that the characteristic time scales correspond to 
break frequencies when calculating the power spectrum 
in terms of w, but simply taking the inverse of the or- 
dinary frequency, 1//, corresponding to the break will 
overestimate tl and tr by a factor of 2tt. 

Physically, a mixed OU process may represent the X- 
ray fluctuations of GBH and AGN under the propagating 
distur bance model, originally developed by iLvubarskil 
(|1997| ). In this model, stochastic disturbances originate 
in the region of the accretion disk outside of the X-ray 
emitting region. These disturbances propagate inwards 
and modulate the X-ray emitting region, which cause 
fluctuations in the X-ray flux. Because many physically 
relevant time scales increase with radius, the X-ray flux 
varies on a broad range of time scales due to the dis- 
turbances originating over a broad range in radii. If 
the disturbances at various radii are independent and 
are well described by an OU process with characteris- 
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tic time scales that increase with radius, then the X-ray 
flux should be well described by a mixture of OU pro- 
cesses. The fact that we can construct a mixed OU pro- 
cess with PSD similar to the X-ray power spectra of GBH 
and AGN is consistent with this interpretation. In this 
case, the break frequencies correspond to the character- 
istic time scales of the disturbances at the maximum and 
minimum radii of the region that contributes to modu- 
lating the X-ray emission. We develop this interpretation 
further in § [3] 

3. MIXED OU PROCESS AS RESULTING FROM 
STOCHASTIC DIFFUSION IN THE ACCRETION DISK 

Similar to us, many authors have attempted to iden- 
tify a particular stochastic process for the accretion rate 
perturbations. While we model the a ccretion rate per- 
turbati ons as following a n OU p rocess . lArevalo fe Uttlev 
2006), iTitarchuk etaLI (|2007l) . and iMisra fe Zdziarski 



2008) have shown that a power spectrum with a 1// 



region can also be created when the accretion rate per 
turbations have a Lorentzian power spectrum, i.e., when 
the disturbances are quasi-periodic. Perturbations of 
this sort are expected if, for example, t he perturbations 
are d ue to Ray leigh - Taylor instability (ITitarchuk et ah 
120071) King et al.l (12004 s ee al so iMaver fe Pringle 
(2006) and IJaniuk fe Czernvl (|2007D 1 model the fluctu- 
ations as resulting from a magnetic field modeled as 
a first order autoregressive p rocess, which is a discret e 
form of the OU process (e.g., iBrockwell fc Davisl 120021) . 
We can build on the work of IWood et al. I (I200lf and 
ITitarchuk et al.l (|2007l ) to show that a mixed OU pro- 
cess can model the viscous response of the accretion disk 
subject to a stochastic driving noise field, for the special 
case when the diffusion equation is linear. 

Conservation of mass and angular momentum result in 
material in the disk obeying a nonlinear diffusion equa- 
tion: 



as 

~dt 



3_9 

r dr 



1/2 



dr 



(r 1 



/ 2 ^e 



(7) 



Here, E is the surface density of the disk and v is the 
kinematic viscosity. Under the assumption that the ro- 
tational frequency of the disk is equ al to the Kepleria n 
value, the accretion rate is given by ( Wood et al. 2001) 



M(r,t) =67rr 1 / 2 ^-(r 1 / 2 z/E). 
or 



(8) 



We can use Equations ([7]) and ((8]) to study the evolu- 
tion of accretion rate perturba tions about the av erage 
steady state solution. Following IWood et al.l ()200lD and 
ITitarchuk et al.1 (|2007f ). define x — r 1 / 2 and a perturba- 
tion to the quantity r x / 2 vYj as u(x,t) — a;A(z/E), and 
therefore AM(r,() cx du(x, t)/dx. The nonlinear diffu- 
sion equatio n for a single perturbation c an now be writ- 
ten as (e.g.. lFrank. King, fc Raine]|2002D 



du 



3 du d 2 u 
AdY.dx 2 ' 



(9) 



When v is independent of E, Equation ^ becomes a 
linear diffusion equation. 

Equation (j9)) only describes the evolution of a single 
perturbation, but in order to create flickering in u(x,t) 



we need to introduce a persistent source of the fluctu- 
ations. We can introduce fluctuations to u(x,t), and 
therefore to AM(r, t), by introducing a stochastic source 
term to Equation (|9|). For simplicity, we only study the 
linear form of Equation ([9]), and therefore assume that 
the viscosity is independent of surface density. While 
it is unlikely that the viscosity is independent of E, it 
is still illuminating to study the linear form of the dif- 
fusion equation in order to better understand how the 
disk responds to a stochastic driving noise, as it permits 
an analytical treatment. Denote the stochastic source as 
8W(x,t)/dt, where dW(x,t)/dt is a spatially dependent 
white noise 0, and therefore W(x, t) is a Wiener random 
field. By definition, dW(x, t)/dt has a flat power spec- 
trum and zero mean. The mean of W(x, t) is also zero, 
and the covariance of W(x, t) is 



E [W(x, t)W(y, a)] = y) mm{t, s), 



(10) 



where E(x) denotes the expected value of x, and <;(x,y) 
defines the spatial covariance of the Wiener random 
field. The Wiener process is essentially a continuous time 
Brownian motion, and therefore the Wiener random field 
W(x, t) may be viewed as an infinite set of Brownian mo- 
tions indexed by the parameter x. If there are no spa- 
tial correlations for W{x,t), then s(x,y) — 5(x — y) and 
dW(x, t)/dt is called a 'space-time' white noise; however, 
in this case, solutions to the linear stochastic diffusion 
equation onl y exist for the case of one spatial dimension 
(|ChowH2007D . 

We assume that the perturbations to the accretion rate 
only exist in a bounded accretion flow defined over a re- 
gion < x < Xmaxi where x max represents the outer 
ed ge of the flow . We adopt the boundary co ndition given 
bv IWood et"aT1 (|2001h and ITitarchuk et alj (|2007f ). with 
u(0,t) — and du(x max ,t)/dt = 0. The first bound- 
ary condition assures that the perturbations of E go to 
zero at the inner boundary of the disk, and the second 
ensures the the disk accretes the perturbations, i.e., the 
perturbations can only leave through the inner bound- 
ary. Including a stochastic source in the linear diffusion 
equation, our stochastic model for the accretion rate per- 
turbations is defined by the following set of Equations: 



du(x,t) 3i/(x) d 2 u{x 7 1) dW(x,t) 



dt 

AM(x,t)=3n 



4 dx 2 
du(x, t) 



dt 



dx 



u(x, 0) = uo(x) 
u(0,i)=0: 



du{a 



dt 



(11) 

(12) 
(13) 
(14) 



Equations (TTT|) -([T4 ^) can be solved using standard meth- 
ods developed for stochastic p artial differen tial equa- 
tions, and a good introduction is iChowl (|2007f ). 

Before presenting the solution, a couple of remarks are 
warranted. First, the model defined by Equations (If II) 
(|I4[) is clearly only an approximation for several reasons. 



2 This notation should not be strictly interpreted as a tradi- 
tional derivative, as the definitions and rules used in stochastic 
calculus are not necessarily the same as those of ordinary calculus. 
However, we ignore this mathematical technicality, as the notation 
8W(x,t)/dt conveys the appropriate interpretation of W(x,t) as 
the integral of dW(x,t)/dt over t. 
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As mentioned before, the actual viscous response of the 
disk will be nonlinear, as the viscosity also likely depends 
on the surface density. In addition, for the noise process 
assumed above, there is nothing to prevent the linear 
model from producing negative values of £ and M, which 
are unphysical. A better approach would be to postulate 
a stochastic model for the viscosity which ensures pos- 
itivity of the solution, and then solve Equation (|9]) nu- 
merically using this perscription for v(r, S, t). However, 
in this work our goal is to study the simpler model, for 
which analytical solutions are possible, in order to better 
understand how the disk responds to a random driving 
noise field, and thus provide some insight on how to con- 
nect features in the power spectra to the physics of the 
accretion disk. Moreover, we will show that the solu- 
tion to the linear diffusion equation is a mixture of OU- 
processes, providing justification for their use in model- 
ing variability. Second, it is important to note that the 
stochastic noise field in the partial differential equation, 
dW(x,t)/dt, is the time derivative of the process which 
is the source of the random perturbations, and is not the 
process itself. And third, the physical interpretation of 
the noise process is left unspecified, however, as many au- 
thors have suggested, dW(x, t)/dt may represent fluctua- 
tions in du{x,t) / dt which a re the result of turbulent and 
chaot ic MHD effects (e.g., lLvubarskiilll997t iKing et all 
2004). Indeed, MHD turb ulence has been observ ed in 
numerical simulations (e.g.. lHawlev fc Krolikl[2001l ) . 

The solution to Equations fTTj) - flT4")) can be expressed 
through the eigenfunctions of the diffusion operator. The 
eigenfunctions, e k (x), and eigenvalues, u) k , are defined by 



Mx) d 2 
4x 2 dx 2 



e k (x) = io k e k (x), 



(15) 



subject to the appropriate boundary conditions. Eigen - 
func tions and eigenva l ues ar e given in lWood et al.l (2001) 
and iTitarchuk et al.l (|2007t ). respectively, for the case 
where the viscosity v(x) is a power-law in x. The eigen- 
functions are orthonormal with respect to the weight- 
ing function p(x) = 4x 2 /(3z/(a;)), and form a complete 
set . The soluti on to Equations (fTT j) - (fT4")) is then given 
bv lChowl ((20071) 



i(x,t) = y\ fc (t)e fc (t), 



(16) 



fc=i 



where u k (t) are a set of stochastic processes. In order to 
determine u k (t), we first note that if the spatial covari- 
ance function of the driving noise, <;(x,y), exists in the 
space spanned by the eigenfunctions, then we can express 
it as an eigenfunction expansion, 



(17) 



fc=l 

where the coefficients are 

4; = / / s{x,y)p{x)e k {x)p(y)e k (y) dx dy. 
Jo Jo 

(18) 

The random field W(x, t) can then be expressed as (e.g., 



ICho^|[2007T) 



oo 

W(x, t) = c, k e k (x)w k (t), 

k=\ 



(19) 



where {w k (t)} are a sequence of independent and identi- 
cally distributed Brownian motions. 

Inserting Equations (fT5"j). p^|) . and (TK)1) into Equa- 
tions (jllj) and (|13p produces an infinite set of ordinary 
stochastic differential equations: 

du k (t) = -co k u k (t)dt + q k dw k (t), k = 1, 2, . . . , (20) 

u k (0)= u a (x)p(x)e k (x) dx. (21) 

Jo 



Comparison of Equations (|2"0"|) and (|21l) with Equation 
([1]) shows that u k (t) are a set of independent OU pro- 
cesses with initial condition u k (0). Therefore, the solu- 
tion to the linear diffusion equation defined by Equations 
CP (HI is 



u(x, t) =]T [ii fe (0)e fc (a;)e- Wfct + s k e k {x)X ou {t, uj k )} (22) 
fe=i 

oc 

AM(x,t)cc^2 



k=l 



u k (0)^Me-^ + ft d -^X ou (t, <m 



dx 



dx 



where X ou (t,uj k ) is an OU process defined by Equation 
(Q} with white driving noise, p, = 0, and <; = 1. The 
first eigenvalue is roughly given by the viscous timescale 
at the outer edge of the disk, T v i S c, depending on how 
v{x) depends on x. In particular, if v{ x] has a power- 
law d ependency on x, then lo\ oc l/r v i sc (jTitarchuk et al.1 
|2007| ). In addition, if the driving noise field is a Gaussian 
process, then the accretion rate perturbations will also 
be a Gaussian process. 

In general, we are interested in the stationary solution, 
which occurs when t — > oo. As a practical matter, this 
occurs after a few viscous times have passed. In this case, 
the initial condition has been forgotten, and the solution 
to the linear stochastic diffusion equation is a mixture of 
OU processes, with the characteristic frequencies given 
by the eigenvalues of the linear diffusion operator, the 
mixing weights given by the product of the eigenfunctions 
of the linear diffusion operator and the coefficients of the 
eigenfunction expansion of the spatial covariance of the 
driving noise field. Then, the covariance function of the 
accretion rate perturbations is 



Cov(AM(x,t),AM{y,s)) = 



Tfc de k (x) de k (y) ^_ L 
j-^ 2uj k dx dy 



•k\t-s\ 



(24) 

and the power spectral density for fluctuations as a func- 
tion of radius is 



4 ( de k {x) 



k=l 



2tt 



dx 



(25) 



The X-ray emission of accreting black holes is thought 
to be released in the inner regions of the disk, and there- 
fore the behavior of the accretion rate fluctuations at 
small x is of primary interest. It is illuminating to in- 
vestigate asymptotic behavior of the power spectrum at 
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x — > for the case where v(x) oc x a . As per the discus- 
sion in § 12.21 the power spectrum of a mixed OU process 
is flat on frequencies shorter than u>x, and therefore the 
fluctuations are decorrelated on time scales longer than 
the viscous time scale at x max . However, on time scales 
short compared to the viscous time scale, the situation 
is more complicated as the power spectrum also depends 
on the spatial covariance structure of the driving noise 
field, ?(x, y), through the set of ft. For simplicity, we first 
assume that there is no spatial correlation in the driving 
noises, and t herefore <?(x, y) = 5(x — y) and ft = 1 for all 
k. Following ITitarchuk et all (j2007h . one can then show 
that the power spectrum of the accretion rate perturba- 
tions at small radii decays as a power-law on timescales 
short compared to the viscous time scale, with the slope 
of the power-law decay steepening when the viscosity in- 
creases less steeply toward larger raddi. Therefore, if 
the viscosity increases more steeply toward higher radii, 
then at small radii the longer time scale perturbations 
are more suppressed relative to the short time scale per- 
turbations. 

In order to investigate the effect that the spatial cor- 
relations in the noise field have on the observed power 
spectra of AM(x = 0, t), we study the case of v(x) oc x 2 . 
In this case the eigenfunctions are si n functions and th e 
eigenvalues are oc (2k - l) 2 fe.g. JWood eiTaTl ^OOF). 
One can show that the power spectrum for AM at small 
radii is then 

P aA (z, w )oc£-£^ (26) 

k=l K 

Wfc = (2fc-l) 2 |. (27) 

As an example, suppose that the spatial correlations for 
W(x,t) decay exponentially in radius with some charac- 
teristic radius r : 

s(x,y) oc e -\ x2 - y2 \ /ra . (28) 

Then, in this case, the coefficients in the eigenfunction 
expansion, c; 2 , are constant for cj^ <C ro/r max and fall off 
as <r| oc 1/w 2 , for oj). S> ro/r max . Inserting this form for 

c; 2 into the Equation for the power spectrum of AM(0, t) 
suppresses the higher frequency modes in the eigenfunc- 
tion expansion, thus supressing the fluctuations on short 
time scales. In particular, for a power spectrum of the 
form Equation (|26[) , a second break will appear near 
uj ro ~ l/T ro , producing a second region which decays 
even faster than the intermediate region. The second 
characteristic time scale, t To , corresponds to the time 
it takes a perturbation traveling at the viscous speed, 
v r ~ r max /T V i SC , to travel across a region corresponding 
to the spatial scale of the driving noise field: 

7"ro ^ TVisc (29) 

This may be the source of the high frequency break 
seen in X-ray power spectra of GBH and AGN, al- 
though it may also be due to an abrupt change in the 
accretion geometry or viscosity (|Churazov et al.1 1200 It 
iPsaltis fc Normanl 120001 ) . Furthermore, we note that if 
the amplitude of the driving noise field increases with ra- 
dius, then ft will decrease with increasing k, steepening 



the slopes in both the intermediate region (uj\ < to < 
ui ro ) and the high-frequency region u> > u> ro . However, a 
decrease in viscosity with increasing radius also has this 
effect, and thus the two situations are degenerate. 

Finally, we note that in order to connect the power- 
spectrum of luminosity fluctations with that of accretion 
rate fluctuations, heat and photon diffusion effects must 
also be considered. The diffusion of heat and photons are 
both governed by diffusion equations, and the radiative 
output of the disk may be described as resulting from 
the convolution of the solution to the heat and photon 
diffus ion equations with th e accretion rate fluctuations 
(e.g.. ITitarchuk et al.l|200 7i), which are stochastic. The 
power spectrum of the radiative output is then the prod- 
uct of the power spectra for the accretion rate perturba- 
tions and the power spectra for the solution to the heat 
and photon diffusion equations. While a detailed treat- 
ment of the effects of heat and pho ton diffusion is beyond 
the scope of this work (but see iSunvaev fc Titarchukl 
ll98CUl985tlTitarchuk et al.ll2007t) . we n ote that based on 
the ab ove discussion and the results of ITitarchuk et all 
(|2007t ). the power spectra of the solution to the linear 
heat and photon diffusion equations will also be of the 
form of Equation (|6|). Therefore, additional breaks will 
occur at frequencies corresponding to the characteristic 
frequencies of the heat and photon diffusion equations, 
and the power spectra of the radiative output will steepen 
even further above these breaks. Thus, the observed fluc- 
tuations in the emission from accreting black holes may 
be interpreted as the result of the response of the disk 
to a random noise field, where the viscous, thermal, and 
radiative response of the disk act as a series of low-pass 
filters on the driving noise field, sequentially suppress- 
ing the variability on time scales short compared to the 
characteristic time scales for viscous, thermal, and pho- 
ton diffusion. Furthermore, if the noise field is spatially 
correlated, then variability on time scales short compared 
to the drift time scale for the characteristic length of the 
spatial correlations is also suppressed. Thus, breaks are 
expected to exist in the power spectra at the frequencies 
corresponding to these time scales, and the slope of the 
power spectra in the different regions contains informa- 
tion on the viscous, thermal, and radiative structure of 
the accretion flow, as well as the structure of the driving 
noise field. 

For completeness, we summarize here the some of the 
results from the above discussion: 

• The mixed OU process is the solution to the lin- 
ear diffusion equation, subject to a random driving 
noise field. The mixing weights are the product of 
the eigenfunctions of the linear diffusion operator 
and the coefficients of the eigenfunction expansion 
of the covariance function of the driving noise field. 
The characteristic frequencies of the OU processes 
are the eigenvalues of the linear diffusion operator, 
subject to the appropriate boundary conditions. 

• The power spectrum of accretion rate perturba- 
tions at small radii, where most of the energy is 
thought to be released, is flat on time scales longer 
than the viscous time scale at the largest radii. On 
timescales shorter than t v ; sc , the power spectrum 
has the form Ci/oj 2 < f AJ y(w) < C2, where C\ 
and C2 are constants. 
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• There is a second break in the power spectrum of 
AM on time scales short compared to the cross- 
ing time of a perturbation traveling at the viscous 
speed across a region of length similar to the char- 
acteristic spatial correlation length of the driving 
noise field. The power spectrum of AM steepens 
on frequencies higher than this break. 

• If the viscosity of the disk increases toward higher 
radii, then in the inner region of the disk the long 
time scale fluctuations are surpressed compared to 
the short time scale fluctuations, thus flattening 
the power spectrum. Similarly, if the amplitude 
of the driving noise increases toward smaller radii, 
the power spectrum also flattens. Thus, these two 
effects are degenerate and cannot be distinguished 
on the basis of power spectra alone. 

• Heat and photon diffusion can introduce additional 
breaks in the power spectra of the emission from ac- 
creting black holes, with the power spectra steep- 
ening above these breaks. 

4. FITTING THE MIXED OU PROCESS 

In the previous section we have suggested an interpre- 
tation of the mixed OU process developed in this work in 
terms of the viscous response of the accretion disk subject 
to a random driving noise field. However, the mixed OU 
process is not limited to linear diffusion models for vis- 
cous evolution of the accretion disk, and in fact one can 
go through the same steps outlined above to show that 
the mixed OU process is a solution for many parabolic 
stochastic partial differential equations, and thus will 
also describe, say, the diffusion of heat in the disk, which 
is governed by the heat equation. However, regardless of 
the physical mechanism, the mixed OU process provides 
an accurate statistical description of the lightcurves of 
GBH and AGN, in the sense that it reproduces much of 
the shape of the power spectra of these objects. As a 
result, estimation of the break frequencies can be done 
directly from the observed lightcurve, and is therefore 
free of the biases due to windowing effects, which are in- 
herent in fitting of power spectra. Some features are not 
captured by this process, such as quasi-periodic oscilla- 
tions and nonlinear effects. An additional quasi-periodic 
component can be modeled as a solution to a set of hy- 
perbolic stochastic partial differential equations, in par- 
ticular the random wave equation. However, inclusion 
of quasi-periodic oscillations to our model is beyond the 
scope of this work. 

4.1. The Likelihood Function 

The mixed OU process may be fit using maximum- 
likelihood or Bayesian techniques. This has the advan- 
tage that the parameters for the process are estimated 
directly from the observed lightcurve utilizing all of the 
information in the data. In particular, the estimates of 
the parameters for the power spectrum model are not 
biased by measurement errors, irregular sampling, or 
other windowing effects due to the finite time span of 
the lightcurve, such as red noise leak. While these bi- 
ases may be corre cted for using Monte Carlo methods 
(jUttlev et al. 2002) , such methods are also computation- 
ally intensive and rely on a \ 2 statistic. Indeed, these 



benefits of time domain modeling are our primary mo- 
tivation in developing the mixed OU process model; we 
sought to obtain a means of estimating the break fre- 
quencies and variability amplitude of X-ray lightcurves 
in a manner that was free of the biases that can affect 
frequency domain techniques. 

In order to calculate the likelihood function of the 
lightcurve, it is necessary to assume a probability dis- 
tribution for the driving noise process. In this work we 
assume that the driving noise dW(t) is Gaussian white 
noise (i.e., W(t) is the Gaussian Wiener process). In the 
context of the linear diffusion equation interpretation, 
this is equivalent to assuming that the driving noise field 
is a Gaussian process. Denote the n observed flux values 
as y\, . . . , un, sampled at time values of tt, ■ ■ • , ijv- The 
observed flux values have independent Gaussian mea- 
surement errors with variances v\ 1 . . . ,vn ■ We denote 
y, t, v to be the vectors containing the N values of flux, 
time, and measurement error variance. In this case it 
is straightforward to calculate the likelihood function, 
as the lightcurves follow a iV-dimensional multivariate 
normal distribution with mean fj, and covariance matrix 
i?y(t) + diag(v), where i?y(t) is given by Equation (0 
calculated at the appropriate values, and diag(v) oper- 
ates on the vector v by constructing a diagonal matrix 
with v along the diagonal. 

Unfortunately, direct calculation of the likelihood be- 
comes computationally prohibitive for most reasonable 
values of N. However, for the discrete mixture the like- 
lihood for an observed lightcurve, modeled according to 
Equation can be efficiently calculated by obtaining 
a state-space representation of the lightcurve, and then 
using the Kalman recursions Q. Using these techniques, 
the likelihood function, p(y|c, u>, fi, <f), is calculated using 
the following recursion formulae: 



p(y|c,w,^,?) = JJ[27rAi] 1/2 exp<j- 

8=1 

-2 



1 (Vi - Vif 

2 A,; 



M 

<r \ - 2 

P3 



3=1 



x, =0 



_2 r 



fli = diag 



cm 



(30) 
(31) 
(32) 

(33) 
(34) 



Aj =diag (exp{-(U - U-i)[ui, . . . , w M ] T ( 



(36) 
(37) 



Xi = AjX^i + "' (yi-i — /J,- c T Xj) 
O i =A < n i _iAf + diag(u oc/!l ) - QtSf/m 

_2 



UOU, 



_ <r 



yi=H + c Xj 
A 2 =c T r2 4 c + Vi 



(39) 

(40) 
(41) 



1 State-space representations and the Kalman recursions are 
commonly used in time ser ies analysis. A good intr oduction to 
these techniques is given by Brockwcll & Davis (2002). 
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Here, c and uj denote the M -dimensional vectors of mix- 
ing weights and characteristic frequencies, respectively, 
x denotes the transpose of x, and Equations (|35|) -(|4T |) 
are only valid for i > 1. Recall that the M character- 
istic time scales are T\ = 1/wi, . . . , Tm = ■ Equa- 
tion (|30[) can then be maximized to calculate maximum- 
likelihood estimates of c, and or used to per- 
form Bayesian inference in combination with a prior dis- 
tribution. Precise determination of the mixing weights 
will in general only be possible with the highest quality 
lightcurves, as the estimates will have a large uncertainty 
with high levels of correlation in their error distribution. 

In practice, an additional constraint needs to be im- 
posed on the norm of c, as it is degenerate with <;. In 
this work we impose a unit squared norm constraint on 
the weights: 

M 

l|c|| 2 =E c ' = l- (42) 

z=l 

Under this normalization, the PSD at frequencies ui 3> 
max{tJi, . . . ,ujm} is 



P(w) 



2tto; 2 ' 



(43) 



from which it is apparent that the amplitude of the high 
frequency PSD is proportional to c; 2 . Thus, ? not only 
gives rate at which variability power is injected into the 
stochastic process, but also defines the amplitude of vari- 
ability at the highest frequencies. The amount of high 
frequency variability in a time interval t% < t < ti is 
given by the integral of Equation (|43[) over this time in- 
terval: 



Var(h <t<t 2 ) = — (i a 

Z7T 



h), t2«l/u H . (44) 



So long as M <C N, calculation of the likelihood via 
Equations (J30j) — ((4TJ) is much faster than direct calcula- 
tion using the autocovariance matrix, Ry(t). In addi- 
tion, we note that Equation ([30| is only valid for Gaus- 
sian measurement errors, and therefore, one must have 
obtained a suitably high number of counts in order to use 
it for X-ray and gamma ray data. An extension to the 
low-count Poisson regime is beyond the scope of our cur- 
rent paper, but will be developed in future work (Kelly 
et al., in preparation). 

Because we are primarily interested in bending power- 
law forms of the power spectrum, we use a mixed OU 
process which can accurately model these forms. We 
have found that the following weighting scheme on a reg- 
ular grid in log uij can provide a good approximation to 
a bending power-law with two breaks: 



M 



-1/2 



_ l-a/2 
■-3 ~ U j 



(45) 



logw, =logw L + 



J-l 

M — I 



(logwtf - logw L ) , j = l,...(M) 



Here, the parameters of the process are n, ?, ujl,ujh, and 
a. The mean of the process is //, <; is the standard devi- 
ation of the driving noise, u>h an d wl are the high and 
low frequency breaks in the power spectrum, and a is 



the slope of the power spectrum in the intermediate re- 
gion. The power spectrum for the process defined by 
Equations (|4"5j) - (l4"6l) is flat for frequencies uj < ujl, de- 
cays as P(u>) oc l/uj a for lol < us < ujh, and decays as 
P(uS) oc 1/ui 2 for uj > ujh. We have experimented with 
different values of M, and find that values of M > 30 are 
more than sufficient; there is little change in the power 
spectrum among values of M > 30. 

The likelihood function may be used to calculate a 
maximum-likelihood estimate. However, in order to 
get reliable estimates of the uncertainties on the mixed 
OU process parameters, we employ a Bayesian approach 
which calculates the posterior probability distribution of 
the parameters, given the observed lightcurve. The prob- 
ability distribution of the parameters given the observed 
lightcurve is 



p(c, uj, n, <r |y) (x p(y |c, w, it, <r)p(c, w, /i, ?) 



(47) 



where the prior distribution is p(c, u, fj., s). For the prior 
under the power-law weighting scheme, we fix the values 
of c and uj to those in Equations (j45|) and (j46|) . and 
assume a uniform prior on — 2 < a < 0, fi, and ? > 0. 
We assume the following prior on the logarithm of the 
break frequencies: 

p(log uj l , log uj h ) = p(log uj h | log ui L )p(\og uj L ) 

1 , N 

= 1 : , UJmin < UJ L < UJ H < W+4S) 

\0gUJ max - log UJ L 

Equation (|48l) is equivalent to placing a uniform prior on 
log ujl over the range ui m i n < ujl < ui max and a uniform 
prior on loguJu over the range ujl < ujh < uj max . The 
upper and lower limits on the characteristic time scales 
that we search for, r max = l/uj min and r mm = l/u; max , 
are chosen to be 10 times the length of the time series, 
and l/10th the smallest time spacing in the time series. 
We use a Markov chain Monte Carlo (MCMC) algorithm 
to obtain random draws of the parameters from Equation 
62]). 

4.2. Assessing the Quality of the Fit 

We suggest two ways of assessing the quality of the 
mixed OU process in fitting the observed lightcurve. 
First, one should analyze the residuals to make sure that 
there are no systematic trends with time, and to ensure 
that there are no significant deviations from normality. 
If there are systematic trends with time, this may imply 
the existence of nonlinear effects. Define the standard- 
ized residuals as 

Xi = (49) 

Under the assumptions used to derive the likelihood func- 
tion, the sequence of Xi should follow a zero mean Gaus- 
sian white noise process with unit variance; i.e., there 
should be no systematic trends in {xt} with time and 
the distribution of the set of Xi should be a Gaussian 
distribution with zero mean and variance equal to one. 
In practice, these assumptions are rarely true, although 
often the Gaussian likelihood approximation is adequate 
for estimating features in the power spectrum, such as 
characteristic time scales and the logarithmic slope of 
the power spectrum, so long as there are not many large 
outliers. 
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The second method for assessing the quality of the fit 
which should be employed is to assess how well the mixed 
OU process fits the observed periodogram. This is simi- 
lar in spirit to the Monte Carlo method of lUttlev et al.1 
(2002), and enables one to assess whether there are fea- 
tures on different time scales which the mixed OU pro- 
cess does not fit well, as such discrepancies may be dif- 
ficult to analyze in the lightcurve residual plot. First, 
one should calculate the periodogram of the observed 
lightcurve. Then, one can simulate a lightcurve hav- 
ing the same time sampling as the observed lightcurve 
for each of the mixed OU process parameters obtained 
from the MCMC sampler. Alternatively, one could also 
just use the maximum-likelihood estimate, but simulat- 
ing lightcurves for each of the MCMC random draws also 
incorporates our uncertainty in the model parameters. 
Then, for each simulated lightcurve one calculates a pe- 
riodogram. These periodograms can then be compared 
with the periodogram for the actual observed lightcurve 
to assess whether there are any time scales where the 
mixed OU stochastic process provides a poor fit. 

5. APPLICATION TO BLACK HOLE LIGHTCURVES 

In order to illustrate our method, and to investigate its 
effectiveness, we applied it to three different data sets. 
The first is an X-ray lighturve of the galactic black hole 
Cygnus X-l in the low/hard state. For this object, the 
data is of high enough quality to clearly see two breaks 
in the power spectrum. The second data set is the sam- 
ple of_l_ocal^AGN_with_l_ong-teTm RXTE data studied 
bv iSobolewska fc P apadakisl (|2009| ) . supplemented with 
XMM lightcurves when available. This sample of ob- 
jects is well-studied, and broken power-law models have 
been estimated and reported by many authors. Both 
of these data sets provide a good test for our statis- 
tical model, as it allows us to compare it with other 
well-established methods. The third data set consists 
of optical lightcurves of AGN from the M ACHO sur- 
vey an d the AGN Watch s ample, taken from lKellv et al.1 
(2009). Kelly et al.l ([20091 ) analyzed these objects using 
only a single OU process, which is not as flexible as the 
mixed OU pr ocess. We exclude the l ightcurves fr om the 
IGiveon et all (|1999Q sample used bv IKellv et all (|2009f) . 
as they are not as well sampled as the MACHO and AGN 
Watch sample. 

5.1. Application to Cygnus X-l 

We used an archival RXTE observation of Cygnus X- 
1 from Oct 23, 199 6 (Obs ID 10241-01-01-000), which 
iNowak et al.l (|1999l) have analyzed and identify as a 
hard-state lightcurve. We extracted the binned PCA 
lightcurve, using all channels, with the standard reduc- 
tion routines. The lightcurve spanned ~ 28 ksec and was 
binned in 0.01 sec intervals. The PSD for this observa- 
tion is shown in Figure [5J where we have subtracted off 
the Poisson noise level. We modeled the PSD using a 
mixture of 32 OU processes, treating the weights as free 
parameters on a fixed logarithmic grid in the character- 
istic frequencies, ujj\ the maximum and minimum values 
of ujj correspond to an order of magnitude shorter and 
longer than the minimum and maximum time spacing 
observed in the lightcurve. 

Unlike for AGNs, the lightcurve has ~ 1.5 x 10 6 data 
points, and computing the likelihood function is ex- 
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Figure 2. The observed PSD for an observation of Cygnus X-l 
in the hard state (data points), compared with the best-fit PSD 
assuming the mixed OU process model (solid line). The mixed 
OU process model provides a good approximation to the observed 
PSD. 

tremely computationally intensive. However, this large 
number of data points has the benefit that the PSD is 
very well determined. As a result, we do not estimate the 
parameters for the mixed OU process from the likelihood 
function of the lightcurve, but instead use least squares 
to fit the model PSD to the observed PSD. The result 
is also shown in Figure O It should be noted that we 
do not bin the lightcurve as an integer number times the 
time-resolution of the instrument mode, and thus some 
aliasing is introduced to the periodogram. However, we 
do not consider this a concern as aliasing is likely to be 
negligible for the low- frequencies used in the fit. As can 
be seen, the mixed OU process is able to provide a good 
approximation to the PSD for this observation of Cygnus 
X-l, although it is not a perfect fit. Further improvement 
could be obtained by making the characteristic frequen- 
cies of the individual OU processes free parameters as 
well; however, we do not do this is this would increase 
the number of free parameters to w 60-70, making the 
least-square optimization computationally difficult. We 
consider this approach of directly fitting the PSD satis- 
factory for Cygnux X-l, even if it is not as powerful as 
maximum-likelihood, as our primary goal in developing 
our method is to estimate features in the PSD, such as 
the break frequencies, from lightcurves that are of signif- 
icantly poorer quality, for which high quality PSDs are 
not available. As such, we apply our model to Cygnus 
X-l merely to illustrate that it provides a good approxi- 
mation of the PSD. 

5.2. AGN X-ray Lightcurves 

5.2.1. Fitting Results 

We applied our method to t he X-ray lightcurves of 10 
local A GN, studied recently bv ISobolewska &: Papadakisl 
(2009) . These objects have high quality lightcurves sam- 
pled over several decades in frequency, and have been 
extensively studied using Fourier-based techniques, mak- 
ing them a good sample for testing our method. The 
sample is summarized in Table [TJ which includes the val- 
ues obtained previously in the literature from Fourier- 
based techniques for the break frequency and logarith- 
mic slope of the power spectra above the break fre- 
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quency. The 2-10 keV RXTE data is presented in 
ISobolewska fc P apadakis (2009), and is important for 
studying the longer time scale behavior. In order to 
also study variations on time scales < 1 day, we included 
XMM lightcurves from the archives, when available. The 
XMM 2-10 keV background subtracted lightcurves were 
extracted using the standard routines of the XMM Sci- 
entific Analysis System, and were binned every 48 sec. 
Both the RXTE and XMM lightcurves were fit simul- 
taneously. The XMM 2-10 keV count rates were con- 
verted to 2-1 keV fluxes using the 2-10 keV flux values 
calculated bv lde Marco et aTl (|2009[ ). with the exception 
of NGC 5506, w here we used the value calculat ed by 
iGuainazzi et all (j2010l) . ISobolewska fc Papadakisl (|2009l ) 
calculated the 2-10 keV luminosities based on the best 
fit spectral model to the 3-20 keV PCA data. 

Our best fit parameters, as well as the 90% confidence 
intervals, are also summarized in Table [TJ In general, 
our values are consistent with values previously obtained 
by other a uthors using t he Fo urier-based Monte Carlo 
method of lUttlev et al.l (|2002[ ). further validating our 
method. However, we obtain consistently smaller errors 
on the estimated parameters. This is likely primarily 
due to the facts that previous work treated the high fre- 
quency logarithmic slope of the power spectrum as a free 
parameter, while in our technique the high frequency end 
of the power law is assumed to decay as l// 2 , and that 
we work directly with the likelihood function. 

In Figure [3] we show the posterior probability distribu- 
tions of the characteristic time scales, logarithmic slope 
of the power spectrum between the break frequencies, 
and the total dispersion in the fractional variations, for 
two sources. The first source, MCG-6-30-15, represents 
some of the best time coverage that we have for our 
sample. The second source, ARK 564, is the only AGN 
with evidence for two break frequencies in its power spec- 
trum, and is thought to be a nalogous to galactic b lack 
holes in the 'very high state' (jMcHardv et all 120071 ): all 
other AGN in our sample are thought to be analogous to 
galactic black holes in the 'high/soft' state. We are able 
to confirm the existence of the two breaks in the power 
spectrum for ARK 564, corresponding to characteristic 
timescales of ~ 300 sec and ~ 2 days. However, for 
MCG-6-30-15, as with all other remaining AGN in our 
sample, there is still considerable probability on values of 
the longer time scale that are longer than the span of the 
lightcurve. As a result, we can only place a lower limit on 
the longer characteristic time scale, as we are only able to 
state with confidence that tl is not less than some value. 
For MCG-6-30-15, we find a lower limit of tl > 37 days. 

We assessed the goodness of fit using the checks de- 
scribed in § 14.21 In general, the mixed OU process pro- 
vided a good fit to the data, although the residuals ex- 
hibited some small deviations from the assumed Gaus- 
sian distribution. We do not consider this a concern, as 
it does not introduce a significant bias in the estimated 
parameters. As an example, in Figure [4] we compare the 
lightcurve of MCG-6-30-15 with the running average for 
the best-fitting mixed OU process. There do not appear 
to be any systematic trends in the residuals accross the 
lightcurve. In addition, in Figure [5] we compare the pe- 
riodogram of the lightcurve for MCG-6-30-15 with that 
expected from the best-fitting mixed OU process, after 
incorporating the uncertainty in the model parameters. 
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Figure 3. Posterior probability distribution for the short and 
long characteristic time scales, PSD slope between the break fre- 
quencies, and dispersion of the X-ray lightcurve for the AGN Akn 
564 (left) and MCG-6-30-15 (right). For MCG-6-30-15 we are only 
able to place a lower limit on the longer characteristic time scale, 
or cquivalcntly, a upper limit on the low-frequency break in the 
PSD. 



We used the Lomb-Scargle periodogram for the RXTE 
data because of its irregular time sampling, while we used 
the standard discrete Fourier transform for the XMM 
data. The mixed OU process, with weights given by the 
bending power-law approximation (Eq. [H]), provides 
a good description of the X-ray lightcurve of MCG-6- 
30-15, being able to model the variability amplitude of 
luminosity fluctuations across a range of time scales. 

As an additional check, we also used a flexible form 
for the mixed OU process, treating the values of the 32 
weights as free parameters, but the break frequencies be- 
ing fixed on a logarithmic grid in uij. This model pro- 
vides more flexible modeling of the PSD, and we fit both 
the Akn 564 and MCG-6-30-15 X-ray lightcurves with it. 
The results are shown in Figure |H1 where we compare the 
'flexible' weighting with the power-law weighting. Note 
that here we have multiplied the PSDs by frequency, for 
easier comparison with previous work. We can still get 
good constraints on the PSD with the flexible weighting, 
although the uncertainties are now larger. For MCG-6- 
30-15, there does not appear to be a significant differ- 
ence in the estimated PSD assuming the flexible weight- 
ing compared to the power-law weighting, although there 
may be some additional 'wiggles' in the PSD which the 
power-law weighting cannot fit. Similarly, the estimated 
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Figure 4. The observed RXTE (left) and XMM (right) 
lightcurves for MCG-6-30-15, compared with the running average 
of the best-fit mi xed OU process model (red line) with power-law 
weig hts (Eq. 1451 ). The standardized residuals, given by Equation 
H49H are also shown in the lower panels. The mixed OU process is 
able to provide a good fit to the X-ray lightcurve for MCG-6-30-15, 
as evidenced by the apparent white noise character of the residuals. 



PSD for Akn 564 assuming the flexible weighting displays 
some wiggles, suggesting that the simple bending power- 
law model is not able to captu r e all o f the features of the 
PSD. Indeed, IMcHardv et all ((20071) find that a sum of 
Lorentzians provides a better fit to the PSD of Akn 564. 
Moreover, 'wiggles' and other deviations from a power- 
law PSD are seen in the X-ray PSDs of GBHs, and thus 
we might also expect them for AGN. 

5.2.2. Comparison with Black Hole Mass 

iSobolewska fc Papadakisl ([20091 ) compiled black hole 
masses for the AGN, most of which a re based on reverber- 
ation mapping dPeterson et al.ll2004D . and we use t he val- 
ues tabulated in ISobolewska fc Papadakisl ([2009). The 
only exception is Mrk 766, which has a mass based on 
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Figure 5. The Periodogram for the X-ray lightcurve of MCG-6- 
30-15 (connected squares) compared with the distribution of peri- 
odograms expected under the mixed OU process model with power- 
law weights (points with error bars). The error bars encompass 
90% of the probability and represent the uncertainty in the ex- 
pected observed periodogram due to uncertainty in the mixed OU 
process parameters, and due to variations in the sampled lightcurve 
due to its stochastic nature. The observed periodogram for MCG- 
6-30-15 is consistent with that expected for the Mixed OU process 
model, confirming that this model can accurately model the vari- 
ability of the X-ray fluctuations of this source across a large range 
in time scales. 



reverberation mapping calculated by Be ntz et~al"1 ([2009) . 
In Figure [7] we plot the characteristic time scale corre- 
sponding to the high frequency break, th, as a func- 
tion of black hole mass. We find a strong correlation 
between the two quantities, although this is not surpris- 
ing as several authors have previously found a correlation 
based on direct PSD fitting of these same sourc es (e.g., 
lUttlev fc McHardvl 120051 IMcHardv etall l2006l) . How- 
ever, using our statistical technique we obtain estimates 
for the break time scales for the two highest mass ob- 
jects, Fairall 9 and NGC 5548, which previous work did 
not include as only lower limits on the break time scales 
were available. This correlation has typically been in- 
terpreted as being driven by a correlation between black 
hole mass and some relevant physical time scale in the 
X-ray emitting region, e.g., the viscous or thermal time 
scale, which corresponds to the suppression of fluctua- 
tions on time scales short compared to this character- 
istic time scale. However, based on the discussion in 
§ [3[ another, possibly related, interpretation is that the 
long time scale break corresponds to the characteristic 
time scale of the physical process driving the fluctua- 
tions, which increases with black hole mass, while the 
correlation between the short time scale break and black 
hole mass is caused by an increase in the characteristic 
spatial scale of the noise process driving the fluctuations 
with increasing black hole mass. 
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Figure 6. PSDs for MCG-6-30-15 (left) and Akn 564 (right) as- 
suming the flexible weighting of the mixed OU process model. Here 
we have plotted fP(f) for better comparison with other work. The 
thin black lines denote 100 random realizations of the PSD from 
its probability distribution, given the observed lightcurve, and the 
thick green lines denotes the median of the random realizations. 
The median value may be considered a 'best-fit' value of the PSD, 
and the spread of the black lines give a graphical representation of 
the uncertainty in the PSD; that probability that the PSD has a 
certain value may be estimated by counting the fraction of black 
lines that intersect that value. The thick red line denotes the best- 
fit PSD assuming the power-law weighting, and the shaded regions 
denote the areas of frequency space not probed by the observed 
lightcurve. While the bending power-law provides a reasonable 
approximation to the PSDs, both sources display additional 'wig- 
gles' in the PSD that are not fit by the bending power-law model, 
similar to what is observed for GBHs. 



We fit a linear relationship between logrn and 
log Man us ing the Bayesian linear regression method of 
iKellvl (|2007l ). which accounts for the measurement errors 
in both log th and log Mbh- It is especially important in 
this case to use a Bayesian method which incorporates 
the measurement errors because we only have 10 data 
points, and thus Bayesian methods are needed in order 
to accurately quantify the uncertainties. We find that on 
average 



t - nn+ 4 - 57 f Mbh 



1.39±0.64 



(ksec) 



(50) 



where the errors are quoted at 90% confidence. We esti- 
mate the intrinsic scatter in th at fixed Mbh to be ~ 0.6 
dex. There is still potentially large intrinsic scatter in th 
at fixed Mbh, which is likely due to scatter in the ac- 
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Figure 7. Dependence of characteristic time scale corresponding 
to the X-ray high frequency break on black hole mass. The error 
bars denote 68% confidence re gion s, and the solid line shows the 
best-fit linear relationship (Eg, 1501 ), Similar to previous work, we 
find a clear dependence of tjj on Mgjj. 



cretion rates at fixed Mbh (e.g.. iMcHardv et alj [2006): 
however, we stress that the estimate of the dispersion in 
the intrinsic scatter is highly uncertain, and more data 
is needed. 

We also searched for correlations of Mbh with the PSD 
slope below the high frequency break a, amplitude of 
driving noise and the total variability amplitude based 
on the mixed OU process model. We did not find any 
significant evidence for trends involving a or the total 
variability amplitude; however, we did find a highly sig- 
nificant anti-correlation between <^ and Mbh, as shown in 
Figured] The trend is surprisingly tight, and implies that 
the fractional X-ray variations on time scales short com- 
pared to the high frequency break are weaker for AGN 
with more massive black holes (see Eq.[44]), or rather 
that the rate at which variability power is injected into 
the lightcurve decreases with increasing Mbh- Within 
the context of the stochastic linear diffusion model, this 
implies that the fractional amplitude of the driving noise 
field, which may be associated with MHD turbulence, de- 
creases with Mbh- The Bayesian linear regression finds 
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-0.79±0.22 



(per cent sec" 



-l/2\ 



(51) 

where the error are quoted at 90% confidence. We esti- 
mate the intrinsic scatter in <; at fixed Mbh to be ~ 0.2 
dex. 

This result is in agreement with previous work which 
has found an anti-correlation between the excess variance 
in the X-r a y ligh t curve and M RH (e.g.. |Lu fc Yul 120011: 
Papadakisl l200l iNikolaiuk et al.1 l2004t lO'Neill et al.1 
2005; iZhou et all 12010(1 . A similar anti -correlation be- 
tween the rate at which variability power is injected 
into the lightcurve and black hole mass has also been 
found for AGN optical variations (KBS09). ICzernv et al.l 
( 2001) find an anti-correlation between Mbh and the fre- 
quency at which the PSD reaches a reference value, which 
is consistent with an anti-correlation between Mbh and 
the amplitude of the high-frequency PSD. However, our 
result differs from most previous work in that we do not 
compute the variance over some range of time scales, di- 
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Figure 8. Dependence of ? on black hole mass; note that ? param- 
eterizes the amplitude of the X-ray PSD on time scales short com- 
pared to the high frequency break, and is proportional to the am- 
plitude of the driving noise field for the mixed OU process model. 
The error bars denote 68% confidence reg ions , and the solid line 
shows the best- fit linear relationship fEq.|50]V There is a clear, 
very tight relationship between ? and Mgjj. 



rectly from the observed lightcurve, but rather we fit a 
parameter q, which is the rate at which variability power 
is injected into the lightcurve. As shown in Equation 
(|44|) . the parameter q 2 is proportional to the short time 
scale variance, but is not strictly equal to it, as it does 
not depend on the time scales probed by the observed 
lightcurve. Indeed, under our adopted normalization for 



the mixed OU process weights (Eq. 



q 2 essentially 



gives the normalization of the PSD on time scales short 
compared to the high frequency break. While the value 
of <; can be precisely determined in an unbiased man- 
ner using our statistical method, this is not the case for 
the observed variance in a lightcurve over some range of 
time scales. The observed variance in a lightcurve over 
some range of time scales is ofte n a poor estimate of 
the actual time-av eraged variance (jVaughan et alJl2003t 
lO'Neill et all l2005). The situation worsens if one is esti- 
mating the variance in the lightcurve using different time 
scales for different sources. 

While we have shown that the ampitude of the high fre- 
quency variance is anti-correlated with Mbh, physically 
interpreting this correlation is difficult. Within the con- 
text of our mixed OU process model the anti-correlation 
implies that the amplitude of the driving noise field, and 
thus the rate at which variability power is injected into 
the lightcurve, decreases with increasing Mbh- However, 
this is also conditional on our adopted normalization for 
the mixing weights (Eq.[42]), and, as we stated earlier, 
the absolute values of the mixing weights are degenerate 
with q 2 . Under our normalization, an anti-correlation 
between <r 2 and Mbh would be expected if the value of 
the PSD in between the break frequencies is indepen- 
dent of Mbh, assuming the PSD slope is equal to —1, 
as u)h is also anti-correlated with Mbh- Thus, the ob- 
served anti-correlation between q 2 and Mbh could be a 
manifestation of the fact that the variability amplitude 
per frequency interval is constant between the break fre- 
quencies, but u>h decreases with increasing Mbh- How- 
ever, if it were true that the ^-Mbh anti-correlation is 
merely an artifact of the th~Mbh correlation, then we 



would expect the <;-Mbh relationship to have a slope of 
—0.5. This is because the parameter q 2 is proportional 
to the value of the PSD at fixed frequency, provided that 
Co 3> ll>h- If the value of u>P(u>) is constant at the high 
frequency break, and if the break frequency scales as 

uh fx -^sif) then we would expect that q oc M 



-1/2 
BH 



However, the dependence of <r on Mbh that we find in 
Equation (51, q oc Mg H 79 ) is steeper than that expected 
under this assumption of a scale-invariant PSD shape. 
This implies that the <;-Mbh anti-correlation is not sim- 
ply an artifact of the th~Mbh correlation, but is real. 

The correlations between the high frequency break or 
? and the black hole mass raise the possibility of using 
the X-ray lightcurve as an additional method for esti- 
mating black hole mass for AGN. It is thus of interest to 
assess the accuracy of the estimates derived from th or 
To do this, we also performed a linear regression of 
log Mbh on \ogrn and log<;, respectively. Note that we 
cannot simply invert Equations (|50l) and (1511) . but need 
to perform the regression after switching the variables. 
Based on our results, mass estimates may be obtained 
as 

\og{M B H/M & ) = (6.75 ± 0.20) + (0.60 ± 0.14) log(rff/[l ksec]) 
\og{M B H/M & ) = (6.36 ± 0.19) - (1.20 ± 0.18) log(?/[(per cent) sec- 
Here, unlike before, we have quoted the errors at 68% 
(l<j) confidence. The dispersion in log Mbh at fixed th 
is ~ 0.4 dex, and the dispersion in log Mbh at fixed g is 
~ 0.2 dex, in agreement with the estimate of iZhou et al.1 
(|20lB when using the total high frequency variability 
amplitude. The probability distribution of the disper- 
sion in mass at fixed th and <j is shown in Figure [9j 
this quantity gives the precision in using these quanti- 
ties as estimates of black hole mass. These results imply 
that mass estimates obtained from the high frequency 
break have similar accuracy to those obtained from the 
broad emission lines (~ 0.4 dex. iVestergaard fc Peterson! 
2006) , but are slightly less accurate than th ose obtained 
from the Mbh-c relationsh ip (~ 0.3 dex, iNovak et al.l 
120061: iGultekin etaLl l2009aT> . The mass estimates ob- 
tained from <;, on the other hand, appear to be the most 
accurate, on average, of all the single-observation tech- 
niques, at least for AGN, and have the advantage that 
the parameter q can be determined from a single X-ray 
lightcurve with significantly higher precision than can the 
break frequency. 

5.3. AGN Optical Lightcurves 

Recent work has suggested that the optical lightcurves 
of AGN are wel l described by a single OU process 
dKellv et al.ll2009t iKozlowski et alJl20iaTMacLeod etafl 
l2010al ). The power spectrum for a single OU process 
becomes flat below the break frequency. It is worth in- 
vestigating whether there is evidence for an intermedi- 
ate region below the break frequency where the power 
spectrum falls off as 1//", (—2 < a < 0), similar to the 
X-ray lightcurves. We can do this by testing whether the 
optical lightcurves are better described by a mixed OU 
process, instead of a single OU process. We used a mixed 
OU process with weights given by Equation (1431) to fit 
the high q uality optical ligh tcurves of 63 AGN from the 
M ACHO (iGeha et all [2003) and AGN Watch samples 
of lKellv et all (|2009l ). The 55 MACHO R- and U-band 
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Figure 9. Probability distribution for the scatter in mass at a 
given th (solid curve) or ? (dashed curve), given the 10 AGN in 
our data set. For comparison, mass estimates obtained from scaling 
relationships involving the broad emission lines and the host galaxy 
properties are ~ 0.4 and ~ 0.3 dex, respectively, implying that c; 
may give the most precise single-observation mass estimates. 

lightcurves typically have seasonal time sampling of ~ 2- 
10 days over ~ 7.5 years. The 8 Seyfert Galaxies from 
the AGN Watch database have optical lightcurves with 
varied time sampling and lengths, and were monitored 
for reverberation map ping. Further deta ils of these sam- 
ples are described in IKellv et al.l ((20090 and references 
therein. 

We searched for objects for which the mixed OU pro- 
cess provided a better fit by finding those objects satis- 
fying the union of the following two criteria: 

• Most of the posterior probability for the ratio of 
break frequencies was found at values u>h /Ve, > 1 ■ 

• Most of the posterior probability for the slope of 
the PSD between the two break frequencies was 
found at values — 2 < a < 0. 

Note that these criteria essentially amount to finding the 
set of PSDs which do not reduce to the Lorentzian shape 
of the single OU process. For most sources, we did not 
find significant evidence that the mixed OU process pro- 
vided a better fit to AGN optical variability than did 
a single OU process. Moreover, the characteristic time 
scales became significantly more uncertain under the 
mixed OU process model, most likely because we allowed 
the slope above the high frequency break to vary. How- 
ever, we did find that for seven of the MACHO sources 
the mixed OU process provided a better fit to both the R- 
and F-band lightcurves. These sources had optical PSDs 
consistent with a single unbroken power-law with slope 
— 2 < a < 0, P(f) oc f a , and thus no break frequencies 
were detected. There were a few additional objects for 
which either the R- or F-band PSD deviated from the 
Lorentzian shape of the OU process, but not both. 

The seven objects for which both the R- and V-band 
lightcurves were better described by the mixed OU pro- 
cess are summarized in Tabled We could not find any- 
thing unusual about these objects other than their timing 
properties. We also fit a flexible form of the mixed OU 
process to some of the better sampled optical lightcurves, 
treating the weights as free parameters on a fixed loga- 
rithmic grid for luj. We could not find any significant 
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Figure 10. Comparison of the optical (solid line, diagonal line 
fill) and X-ray (dashed line, solid fill) PSDs for the four objects 
in our sample for which we have both good X-ray and optical 
lightcurves. The lines denote the best-fit PSDs and the shaded 
regions denote the 68% confidence regions for the PSDs, assuming 
the mixed OU process model with power-law weighting. In general, 
the optical PSDs fall off faster toward higher frequencies, implying 
that the short time scale variations are more suppressed in the op- 
tical lightcurves. However, the optical lightcurves appear to have 
more variability on longer time scale than the X-ray lightcurves, 
implying that the optical variations are not solely due to reprocess- 
ing of the X-ray emission. The only exception is Akn 564, which 
has a high accretion rate and is thought to be analogous to galactic 
black holes in the very high state. 



evidence that the optical PSDs diverged from either a 
single power-law or Lorentzian PSD. 

We have both optical and X-ray lightcurves for Akn 
564, NGC 3783, NGC 4051, and NGC 5548, and it 
is worth comparing the two estimated PSDs for each 
source. In Figure [10] we compare the estimated opti- 
cal and X-ray PSDs for each source assuming the mixed 
OU model with power law weighting. For all sources we 
subtracted off the host galaxy flux, w here we used the 
values reported by iBentz et al.l 1)20090 for NGC 3783, 
NGC 4051, and NGC 5548, and the value reported by 
iShemmer et al.l (|2001l ) for ARK 564. In general, the op- 
tical fluctuations are less variable then the X-ray fluctu- 
ations, except possibly on time scales > 100 days. Fur- 
thermore, the optical PSDs fall off steeper toward high 
frequencies than the X-ray PSDs, and thus the short time 
scale optical fluctuations are more strongly suppressed 
relative to the X-ray fluctuations. These results are con- 
sistent with previous comparisons of optical and X-ray 
PSDs (e.g., ICzernv et al.)ll99l 120031: lUttlev et al.[ l200a 



lArevalo et af.ll2008l 120091: IBreedt et afl l2009. 2013). 

One possible interpretation of these results is that the 
different PSDs represent how the different media in the 
two regions respond to a common driving process that 
is turbulent in both space and time, with the region 
emitting the optical photons more strongly suppressing 
the short-time scale noise fluctuations. This is consistent 
with the fact that ? is anti-correlated with Mbh for both 
the X-ray and optical lightcurves. For example, within 
the propagation class of models, this implies that while 
the source of the accretion rate perturbations may be 
the same, the optical emitting region more strongly sup- 
presses the short time scale accretion rate perturbations, 
but when these perturbations reach the X-ray emitting 
region (i.e., the corona), which may also be subject to a 
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driving noise field, the short time scale fluctuations are 
not as strongly filtered as this medium is not as 'stiff'. 
However, if the optical PSDs continue to rise above the 
X-ray PSDs at longer time scales, as it appears they 
do, then the total amplitude of variability in the opti- 
cal fluctuations surpasses that in the X-ray fluctuations. 
This is a departure from the behavior seen in soft state 
GBH, where the disk emission is significantly less vari- 
able than the corona, suggesting a descrepancy between 
the spectral-timing properties of GBHs and AGN. In ad- 
dition, larger optical variability amplitude than X-ray 
variability amplitude implies that either the source of 
the variations is different for the two emission regions, or 
that the power in the disturbances is somehow dissipated 
such that their variability amplitude has decreased by the 
time they reach the X-ray emitting region. However, the 
results also imply that all of the optical variability can- 
not be due to reprocessing of X -ray emission, as has been 
shown in previous work (e.g., iCzernv et afl[T99l [20031: 
lUttlev et al.ll2003i ). The only possible exception appears 
to be Akn 564, where the X-ray variability amplitude is 
larger at all time time scales. Akn 564 is noteworthy as it 
is the only AGN known to exhibit a second low-frequency 
break, and the additional X-ray variability may be due 
to an increase in the strength of the X-ray power-law 
component, such as is seen in GBH in the VHS. 

6. SUMMARY 

In this work we have developed a new model for the 
lightcurves of accreting black holes. The model is based 
on a mixture of Ornstein-Uhlenbeck stochastic processes 
and also happens to be the solution to the linear diffu- 
sion equation perturbed by a spatially correlated noise 
field. In essence, our model for quasar lightcurves may 
be thought of as a type of basis expansion, where the 
basis for the lightcurve is a set of independent stochastic 
processes. This allows flexible modeling of the PSD so 
long as — 2 < d log P/d log / < for all frequencies of in- 
terest, and we show that the mixed OU process provides 
a good approximation to the X-ray lightcurves of GBHs, 
and the optical and X-ray lightcurves of AGN. Our main 
results are 

• The mixed OU process is the solution to the linear 
diffusion equation perturbed by a spatially corre- 
lated noise field, and thus the mixed OU process 
describes the evolution of viscous, thermal, or ra- 
diative perturbutions due to a driving noise, to the 
extant that the viscous, thermal, or radiative re- 
sponse of the accretion flow follows a linear diffu- 
sion equation. Within this interpretation, the low 
frequency break in the PSD corresponds to the dif- 
fusion time scale in the outer region of the accretion 
flow, and the shape of the PSD above the low fre- 
quency break depends on the viscosity of the flow 
(or, more generally, the diffusion coefficient). If the 
noise field is spatially correlated, then an additional 
high frequency break in the PSD may exist at the 
frequency corresponding to the time it takes a per- 
turbation traveling at the viscous speed to cross 
the characteristic spatial scale of the noise field. 

• We derive the likelihood function of a lightcurve 
for the mixed OU process, given by Equations (l30l) - 
(l4Tj) . This allows one to estimate the parameters of 



the mixed OU process for an observed lightcurve, 
and equivalently the PSD, in a manner which is 
computationally efficient, unbiased by red noise 
leak and aliasing, and fully accounts for irregular 
sampling and measurement errors. If the parame- 
ters for a bending power-law are desired, then we 
also derive a form for the mixing weights which ap- 
proximates a bending power-law, given by Equa- 
tion fljgj) . 

• We show that the mixture of Lorentzians PSDs im- 
plied by the mixed OU process stochastic model is 
a good fit to the PSD of an observation of Cygnus 
X-l in the hard state. 

• We applied our mixed OU process model to the X- 
ray lightcurves of 10 well-studied local AGN and 
show that it provides a good description of their 
lightcurves, and we recover many of the PSD fea- 
tures which have been obtained previously from di- 
rect fitting of the observed PSDs. We use the esti- 
mated break frequencies to recover the correlation 
between time scale of the high frequency break and 
black hole mass seen previously with this sample; 
this includes new estimates of the high frequency 
time scales for Fairall 9 and NGC 5548, for which 
previously only upper limits were obtained. 

• We find a significant anti-correlation between black 
hole mass and the amplitude of the driving noise 
fluctuations, which also parameterizes the ampli- 
tude of the high frequency PSD. The form of this 
correlation is similar to what has been observed 
in optical lightcurves (KBS09), suggesting that the 
origin of the optical and X-ray variability for AGN 
is at least partially shared. This correlation is even 
tighter than that between Mbh and the break fre- 
quency. Mass estimates obtained using ^ appear to 
have an uncertainty of ~ 0.2 dex, potentially mak- 
ing it the best single-observation mass estimator 
for AGN. 

• We did not find significant evidence that most opti- 
cal lightcurves of AGN deviated significantly from a 
single OU process, i.e., a Lorentzian PSD. However, 
we did find evidence that the optical lightcurves of 
7 AGN in a sample of 55 AGN from the MACHO 
survey exhibited PSDs which were flatter than the 
l// 2 shape seen in the majority of the sample, 
but the origin of this difference is unclear. These 
sources may exhibit flatter PSDs because we are 
seeing them on the 1// Q part of the PSD, i.e., the 
break time scales are shorter and longer than the 
minimum and maximum time scales probed by the 
lightcurves. 

• We compare the optical and X-ray PSDs for Akn 
564, NGC 3783, NGC 4051, and NGC 5548 and 
find that for the three NGC sources the short time 
scale fluctuations are more suppressed in the op- 
tical, but the long time scale fluctuations appear 
to be stronger in the optical. This implies that 
the total variability amplitude in the optical fluc- 
tuations is larger than that in the X-ray fluctua- 
tions, and thus the optical fluctuations cannot be 
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caused solely by reprocessing of X-ray emission. In 
addition, if both the optical and X-ray fiuctations 
are caused by inwardly propagating fluctuations in, 
say, the accretion rate, then some of the power in 
the fluctuations must be dissipated by the time 
they reach the X-ray emitting region. Akn 564, 
the only AGN believed to be in the very high ac- 
cretion rate state, does not exhibit this behavior, 
and the optical fluctuations are always less variable 
then the X-ray fluctuations, at least over the range 
of time scales we probe. 
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APPENDIX 
LEVY PROCESSES 

Levy processes form a general set of processes that include Brownian motion and compound Poisson processes, 
both relevant in astrophysics. Basically, a stochastic process is a Levy process if it has stationary and independent 
increments. Brownian motion, or rather the Wiener process, is a random walk process with a power spectrum that 
decays as 1/ui 2 . A white noise process is the derivative of Brownian motion, and is a stochastic process with mean 
zero and a flat power spectrum; conversely, Brownian motion is the integral of white noise. A Brownian motion for the 
process driving the luminosity fluctuations might be considered a reasonable approximation if the input driving noise, 
dW(t), has a decorrelation time scale short compared to the characteristic time scale of the luminosity fluctuations, 
therefore implying that dW(t) resembles white noise. Alternatively, a compound Poisson process may be considered to 
describe the variability in luminosity caused by flares above an accretion disk. If the number of flares in a time interval 
follows a Poisson distribution, and the luminosities of the flares are randomly drawn from some other probability 
distribution, then the sum of the flare luminosities is a compound Poisson process. In general, most of the results from 
§[2] are valid for any Levy process with zero mean and unit variance. The assumption that the driving noise, dW(t), 
has zero mean and unit variance is not very restrictive, as [i and s can always be rescaled to make this so. 

CONTINUOUS MIXTURES OF ORNSTEIN-UHLENBECK PROCESSES 
A continuous mixture of independent OU processes, denoted as Y(t), may be obtained as 

/•oo />oo 

Y (t)=ft+ / X(t,wo,OC(wo,0 du> tfc (Bl) 
Jo Jo 

Here, C(wo,<r) is a weighting function which describes how much a given OU process contributes to the mixture, and 
the notation X(t,cvo,<;) defines an infinite set of independent OU processes indexed by the value of ojq and we 
will refer to C(u>o,<;) as the 'mixing function'. We assume that the individual OU processes have zero mean, and 
that the mean of the mixed OU process is /i. Some previous work has focuse d on the special case where C(cjo,<t) 
is a probability distributio n (e.g.. Ilgloi fc Te rdik 1999; Eliazar & Klafter 2009). In general, the process Y(t) is not 
Markovian but stationary (jEliazar & Klaftcr H2009f K 

We are allowed considerable flexibility in modeling the lightcurve Y(t) via the mixing function, and as an example 
we focus on the special case where C(u;o,<r) is a step function over some range of ojq. Specifically, 

Ci - UJ °' ^ = {k otherwise" ^ ' ( B2 ) 

The aut ocovariance function and power spectrum can be calculated for this choice of mixing function using the formulae 
given in lEliazar fc Klafterl (|2009t ). The autocovariance function is 

R Y (t) = ^{E 1 (t/T L )-E 1 (t/T H )}, (B3) 

where tl, = 1/ujl and th — ^/ojr are the maximum and minimum characteristic time scales of the mixture, and E\{x) 
is an exponential integral of first order: 

El (x) = [°° exp t-^ d y, x>0. (B4) 
Ji V 
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The variance of the continuous mixed OU process, Y(t), is calculated by taking the limit of Equation (|B3[) as t — » 0: 

Var[Y(t)} = lim R Y (t>) = q - In f ^) . (B5) 

Because the individual OU processes involved in the mixture are assumed to have zero mean, the mean of Y (t) is fi. 
The power spectrum for the continuous mixed OU process with weighting function given by Equation (|B2|) is 



Three regions are of interest , 
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(B6) 



(B7) 



Therefore, the mixed OU process resembles white noise on time scales St ^ tl, 'pink' or 'flicker' noise on time scales 
Th <C St <C Tl, and red noise on time scale St <C Tff. For a lightcurve that follows a mixed OU process, tx represents 
the maximum time scale on which Y(t) is correlated; on time scales longer than tx the lightcurve 'forgets' about its 
previous behavior. 

We can also formulate a mixed OU process with common driving noise. This process was studied bv lEliazar & Klafter 
(2009). Physically, we might expect this type of process to arise if, for example, the fluctuations in the optical lightcurve 
were the result of reprocessing due to a common source which irradiates a range of radii, with the reprocessing time 
scales increasing with increasing radii. Denote this type of mixed OU process as Y(t). In this case, the power spectrum 
using the mixing function defined by Equation (|F32|) is 
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The mixed OU process with common driving noise exhibits less long-time scale variations than does the mixed OU 
process with independent driving noise, due to the correlation in the processes created by a common driving noise. 
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J The values reported by McHardy et al. 12004) are for a sharply broken power-law, as the authors did not state confidence regions for the bending power-law PSD model. 



Table 2 

MACHO Sources with Optical Lightcurves that Deviate from a Lorentzian PSD 
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